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A discrete multistate kinetic model for water-wire proton transport is constructed and analyzed 
using Monte-Carlo simulations. The model allows for each water molecule to be in one of three states: 
oxygen lone pairs pointing leftward, pointing rightward, or protonated (HsO + ). Specific rules for 
transitions among these states are defined as protons hop across successive water oxygens. We then 
extend the model to include water-channel interactions that preferentially align the water dipoles, 
nearest-neighbor dipolar coupling interactions, and coulombic repulsion. Extensive Monte-Carlo 
simulations were performed and the observed qualitative physical behaviors discussed. We find the 
parameters that allow the model to exhibit superlinear and sublinear current-voltage relationships 
and show why alignment fields, whether generated by interactions with the pore interior or by 
membrane potentials always decrease the proton current. The simulations also reveal a "lubrication" 
mechanism that suppresses water dipole interactions when the channel is multiply occupied by 
protons. This effect can account for an observed sublinear-to-superlinear transition in the current- 
voltage relationship. 
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INTRODUCTION 

The transport of protons in aqueous media and across 
membranes is a fundamental process in chemical re- 
actions, sol vation, and pH reg u lation in cellular en - 
vironments [Alberts et al 1994 iGrabe fe Oster 200l| . 
Proton transport in confi ned geome tries is also rele- 
vant for ATP synthesis iBover 1 997| and light trans- 
duction by bacteriorhodopsin jLanx^W9^. In this 
paper, we develop a lattice model for describ- 
ing proton transport in one- dimensional environ- 
ments. This study is motivated by numerous mea- 
surements of proton conducti on across channels em- 
bedded in lipid membranes | Akeso n fc Deamer 1991 



Busath et al. 1998llCotten et al 19991 ICukierman 1997 



Dea^ieTT^H^J^smim^Tr^I'T^^^ Experiments are 
typically performed using membrane-spanning grami- 
cidin channels that are only a few Angstroms in diame- 
ter. This geometric constraint imposes a single-file struc- 
ture on the configurations of the interior water molecules 
|Hille fc Schwarz 19781 Iffladkv fc Havdon 1972| . 

Under the same electrochemical potential gradients, 
conduction of protons across ion channels occurs at a 
rate typically an order of magnitude higher than that of 
other small ions. This supports a "water- wire" mech- 
anism lAkeson fc Deamer 199ll iNaele fc Horowitz 1975 
iNagle fc Tristan-Nade 19831 iNagle 19871. fir st proposed 
by Grotthuss [Agmon 19951 iGrotthuss 18 06], Across a 
water- wire, protons are shuttled across lone pairs of water 
oxygens as they successively protonate the waters along 
the single-file chain. However, since the hydrogens are 
indistinguishable, any one of the hydrogens in a water 
cluster (e.g., any of the three hydrogens on a hydronium) 
can hop forward along the chain to protonate the next 
water molecule or cluster of water molecules (cf. Fig. [Q. 
This mechanism naturally allows much faster overall con- 



duction of protons compared to other small ions which 
have to wait for the entire chain of water molecules ahead 
of it to fluctuate across the pore in order to traverse the 
channel. 

A peculiar feature of measured current-voltage re- 
lationships is a crossover from sublinear to super- 
linear behavior as the pH of the reservoirs is low- 
ered. Measurements by Eisenman [Eisenman et al. 1980j 
were carried out in symmetric solutions in the 1-3 pH 
range, and the results were rece ntly reproduced by 
Busa th et al. [Busatli ^ al 1998j and Rokitskaya et 
al [Rokitskava et al. . 2 002J . These experiments were 
performed using simple, relatively featureless gramicidin 
A (gA) channels. One leading hypothesis is that the 
nonlinear proton current-voltage relationships arise from 
the intrinsic proton dynamics within such simple chan- 
nels. Specifically, multiple proton occupancy and re- 
pulsion among protons within the channel may give 
rise to the observed nonlinearitv | H ille_fc_ Schwarz 19781 
Phillips et al. 19991 ISchumaker et al. 200lj . 

There have been a number of recent theoretical 
studies of water-wire proton conduction. Extensive sim- 
ulations on the quantum dynamics of proton exchange 
in essentially small, representative water clusters in 
vacuum have been used to predict microscopic hop- 
ping rates between w ater clusters iBala et al. 1994 , 
Sadgeghi fc Cheng 19991 iMarx et al. 1999 . 

Mavri fc Berendsen 1999, iMei et al. 1998 , 

Sagnella et al. 19961 ISchmitt fc Voth 1999j . Pomes and 
Roux |p'oTn^^^Rpux^^96| have performed classical 
molecular dynamics (MD) simulations on water-channel 
interactions, proton hopping, and water reorientation. 
They derive effective potentials of mean force describing 
the energy barriers encountered by a single proton within 
the pore. Since MD simulations are presently limited to 
only processes that occur over a few nanoseconds, none 
of these computational methods are efficient at probing 
very long time, steady-state transport behavior. On 
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a more macroscopic, phenome nological level, Sagnella 
et al. ISaenella fc Voth 19961 and Schumaker et al. 
[Schumaker et al 200(1 ISchumaker et al 200 lj have 
considered a the long-time behavior of a single proton 
and dipole "defect" diffusing in a single- file channel. 
The parameters used in these studies, including effective 
energy profiles and kinetic rates, were derived from MD 
simulations. Although the basic underlying structure 
assumed by all of these transport models qualitatively 
resembles the Grotthuss mechanism, they have not 
addressed multiple proton occupancy. 

In this paper, we will explore the intrinsically non- 
linear proton dynamics along a single-file water-wire. 
We use a dynamical lattice model that defines the dis- 
crete structural states of the water-wire that approxi- 
mate the continuous molecular orientations. Although 
the lattice model provides a different approach from MD 
simulations, it is more amenable to analysis at longer 
time scales, yet is connected to the microphysics inher- 
ent in MD simulations provided a consistent correspon- 
dence between the parameters is made. Rather than 
enumerating all possible molecular configurations, our 
lattice appr oach is resembles that deve loped for molecu- 
lar motors |Fisher fc Kolomeiskv 1999|. mRNA transla- 
tion IMacDonald fc Gibbs 19691 IChou 20031. traffic flow 
|KarimipoTu^i999l ISchrecke^u3en^^t^l99a . and ion 
and water transport in single-fi le channels [Chou 19981 
IChou 19991 IChou fc Lohse 1999| . Here, the proto n occu- 
pancy along the water- wire will be self-consistently deter- 
mined by the prescribed lattice dynamics. The parame- 
ters used in our model are transition rates among discrete 
states that in principle can be independently computed 
from relatively short-time MD simulations. Despite the 
approximations inherent in our discrete model, it quali- 
tatively treats the effects of proton-proton repulsion and 
water-water dipole interactions. 

MODEL AND METHODS 

Qualitatively, protons hop from oxygen to oxygen dur- 
ing transport. The successive hops clearly do not have 
to involve an individually tagged proton; in this respect, 
proton currents resemble electrical conduction in a con- 
ductor. Many measurements of proton conduction across 
membranes are performed on the gramicidin model sys- 
tem. The interior diameter of gramicidin A (gA) is ~ 3—4 
Angstroms and can only accommodate water in a single 
file chain. Although the number of water molecules in 
this chain is a fluctuating quantity, their dynamics in and 
out of the channel will be assumed to be much slower than 
that of t heir orientational rearrangements and proton 
hopping [Hummer et al. 200ll iKalra et al. 2003| . We 
can thus treat the water wire as containing a fixed, av- 
erage number of water molecules. Within typical trans- 
membrane channe l s, are 7V 8 — 2 6 single-file waters 
[Levitt etal. 197SLIWu fc Voth 2003| . 

Figure 1A shows a schematic of our model. We 
first assume that each "site" along the pore is occu- 
pied by a single oxygen atom which may either be 



part of neutral water, H2O, or a hydronium (H30 + ) 
ion. Although protonated oxygens in bulk are often 
associated with larger complexes such as HsO^~ (Zun- 
del cation), or HgO^ (Eigcn cation), in confined ge- 
ometrie s, the formation of the larger complexes is sup- 
pressed [Lvnden-Bell fc Rasaiah 1996], Furthermore, we 
will show that our discrete model depicted in Fig. 1A 
can incorporate the dynamics of reactive proton transfer 
among transient clusters by an appropriate redefinition 
of a lattice site to contain the entire cluster. 
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FIG. 1: (A) Schematic of an N = 11, three-species exclusion 
model that captures the steps in a Grotthuss mechanism of 
proton transport along a water wire. For typical ion channels 
that span lipid membranes, N ~ 10 — 20. The transition rates 
are labeled in (B) and in the legend. Water dipole kinks are 
denoted by thick lines. 



Neutral waters have permanent dipole moments and 
electron lone-pair orientations that can rotate thermally. 
For simplicity, we bin all water dipoles (hydrogens) that 
point towards the right as "+" particles, while those 
pointing more or less to the left are denoted "— " par- 
ticles. The singly protonated species H30 + is hybridized 
to a nearly planar molecule. Therefore, we will assume 
that hydronium ions are symmetric with respect to trans- 
ferring a proton forward or backwards, provided the ad- 
jacent waters are in the proper orientation and there are 
no external driving forces (electric fields). Each lattice 
site can exist in only one of three states: 0, +, or — , corre- 
sponding to protonated, right, or left states, respectively. 
Labeling the occupancy configurations <ii — { — 1, 0, +1}, 
allows for fast integer computation in simulations. 
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In addition to proton exclusion, the transition rules are 
constrained by the orientation of the waters at each site 
and are defined in Fig. IB. A proton can enter the first 
site (i = 1) from the left reservoir and protonate the first 
water molecule with rate a only if the hydrogens of the 
first water are pointing to the right (such that its lone- 
pair electrons are left-pointing, ready to accept a proton) . 
Conversely, if a proton exits from the first site back into 
the left reservoir (with rate 7), it leaves the remaining 
hydrogens right-pointing. In the pore interiors, a pro- 
ton at site i can hop to the right (left) with rate p+(p_) 
only if the adjacent particle is a right(lcft)-pointing, un- 
protonated water molecule. If such a transition is made, 
the water molecule left at site i will be left (right) point- 
ing. Physically, as a proton moves to the right, it leaves 
a wake of — particles to its left. A left moving proton 
leaves a trail of + particles to its right. These trails of — 
or + particles are unable to accept another proton from 
the same direction. Protons can follow each other suc- 
cessively only if water molecules can reorient such that 
these trails of +'s or — 's are thermally washed out. Wa- 
ter reorientation rates are denoted k± (cf. Figs. IB and 
2). Protons at the rightmost end of the water wire (at 
site i = N), exit with rate /3, which is different from p + 
since the local microenvironment (e.g., typical distance 
to acceptor electrons) of the bulk waters that accept this 
last proton is different from that in the pore interior. 
From the right reservoir, protons can hop back into the 
water wire with rate S if a water in the "— " configura- 
tion is at site i = N. The entrance rates a and 6 are 
functions of at least the proton concentration in the re- 
spective reservoirs. Figure 2 shows a representative time 
series of the evolution of a specific configuration. The 
rate-limiting steps in steady-state proton transfer across 
biological water ch annels are thought to be associated 
with water flipping [Pomes fc Roux 1998| . 

The lattice discretization for individual H3O 4 " ions 
need not be interpreted literally. Larger complexes can 
be effectively modeled by reinterpreting p± , k± , and the 
basic unit of hopping for the proton. For example, if 
certain conditions obtain, where ions are predominantly 
two-oxygen clusters (H5OJ), we defined each pair of wa- 
ters as occupying a single lattice site, k± as an effec- 
tive reorientation time for the following pair of waters, 
and p± as the hopping rate to an adjacent oxygen lone- 
pair. The Grotthuss water-wire mechanism is qualita- 
tively preserved as long as the proper identification with 
the microphysics is made. 

All eight "parameters" used in our model (the 
rates p±, k±, a, /?, 7, 5), can be related to measured 
bulk quantities or derived from short-time MD simu- 
lations. They are a minimal set and are equivalent 
to the numerous bulk par ameters used in other mod- 
els [Schumaker el al. 2001) . such as the bulk proton dif- 
fusion constant, water orientational diffusion constants, 
etc. Using similar MD approaches then, one should be 
able to approximately fix the parameters used in our 
model. For example, variations in the potential of mean 
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FIG. 2: A time series depicting a number of representative 
transitions obeying the dynamical constraints of our model. 
A proton (0) at site i can move to the right with rate p+ 
only if site i + 1 is occupied by a properly aligned (lone-pair 
electrons pointing to the left) water molecule (+). When 
a proton leaves site i to the right, it leaves behind a water 
in state "— ", with lone pair electrons pointing to the right. 
Protons at site i can also move to the left with rate p_ if 
site i — 1 is a water in the "— " state. In this case, a water 
is left behind at behind site i in the "+" state. The neutral 
water molecules must flip (+ «-> — ) in order for a nonzero 
steady-state current to exist. 
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FIG. 3: (A — D) Energy differences between final and initial 
states which involve a change in ferroelectric coupling, net 
dipole moment, and repulsive interactions. (E) A representa- 
tive energy barrier profile for H = K — R = V = (dashed 
curves). The energy profile for H, K, R, V 7^ for a transition 
between the states considered in (D) is shown by the thick 
solid curve. 
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force along the pore (resulting from interactions of the 
different species with the constituents of the pore in- 
terior) are embodied by site-dependent transition rates 
p± and k±. Thus, MD-derived potentials of mean force 
used in previous models can also be implemented within 
our lattice framework. Such effects of local inhomo- 
geneities in the hopping rates have been studied an- 
alytically and wit h MC simulations in related models 
|Kolomeiskv 1998| . 

The basic model described above has been studied an- 
alytically in certain limits where exact asymptotic re- 
sults for the ste ady-state proton current J were de- 
rived |Chou 2002j . However, this study did not explic- 
itly include any interactions other than proton exclusion 
and proton transfer onto properly aligned water dipoles. 
Effects arising from forces such as repulsion between 
protons in close proximity, interactions between water 
dipoles and external electric fields, and dipolar coupling 
between neighboring waters need to be considered. 

In Fig|3j4, a proton moves down the electric potential 
reducing the total enthalpy by V, and a right-pointing 
dipole is converted into a left-pointing dipole at an en- 
ergy cost of H . Since both initial and final states have 
adjacent, repelling protons, the repulsion energy R does 
not enter in the overall energy change. In Fig. |3|B, a 
proton moves down the potential (— V), a "+" water is 
converted to a "— ", (+H), a dipole "domain wall" is re- 
moved (-K), and the repulsive energy between adjacent 
charged protons is relieved (— R). The representation of 
these nearest neighbor effects can be succinctly written 
in terms of the energy of a specific configuration 

N-l N 
E[{ai}] =-K^ Wi+i -HY,^ + 

i=l i=l 

«EiIi X (l - *?)(! - °f+i) - V££x i(l - 

(l) 

The H,K,R,V parameters used in are all in 

units of fcsT and represent 

• H: energy cost for orienting a water dipole against 
external field 

• K: energy cost for two oppositely oriented, adja- 
cent dipoles 

• R: repulsive Coulombic energy of two adjacent pro- 
tons 

• V: energy for moving a charged proton one lattice 
site against an external field. 

V as the change in potential that a proton incurs as it 
moves between adjacent waters. The total transmem- 
brane potential V mem brane = NV. The local dielectric 
environment across a channel ca n induce a spatially 
varying effec tive potential Vi <^<jv |Edward^^z^^002 , 

Jordan 1984 IPartenskii fc Jordan 1992 , 

Sveanow fc von Kitzing 1999| . In this study, we 



neglect this variation and assume constant V across the 
lattice. 

In order to connect the quantities H, K, i?,and V to 
the rates a, 0, 7, 5,p±, k±, we will assume the transitions 
occur over thermal barriers. Although barriers to pro- 
ton hopping may be small, we employ the Arrhenius 
forms in order to obtain a simple relationship so that 
qualitative aspects of the effects of H, K, R, and V can 
be illustrated. Activation-energy-based treatments for 
conduction acr oss gramicidin channels have be en pre- 
viously studied [Chernvshev fc Cukierman 2002J. When 
the more complicated interactions and external poten- 
tials are turned on, the effective transition rates £ = 
{a, (3, 7, S,p±,k±} on which we base our Metropolis 
Monte-Carlo become 



e = £oexp(_j, (2) 

where £0 = {&o, Po, Jo, So,Po, ko} are rate prefactors 
when H, K, R, V and AE are zero. In defining Eq. |3 
we have assumed that the energy barrier due to the 
difference AE = E^}] - E[{(7i}} {{a[} and {ctJ are 
the final and initial state configurations, respectively) is 
evenly split between the barrier energies in the forward 
and backward directions. We use the convention that 
p + = P- = po and k + = k- = ko when V — and H = 0, 
respectively. The constraints and the state-dependent 
transition rates determined by Eqs. 1 and 2 completely 
define a nonequilibrium dynamical model which we study 
using MC simulations. Note that in the original model 
(Fig. 2) we do not assume transition barriers, but rather 
only that the dynamics are Markovian. 

We first gained insight into the dynamics by consider- 
ing numerical solutions to the full Master equation for a 
short three site (N = 3) channel. If we explicitly enu- 
merated all 27 = 3 3 states of the three site model, the 
Master equation for the 27 component state vector P is 

^f = MF(t), (3) 

where M is the transition matrix constructed from the 
rates £. In steady-state, the P, are solved by inverting 
M with the constraint X^ =1 Pi = 1. The steady-state 
currents are found from the appropriate elements in Pj 
times the proper rate constants in the model. For exam- 
ple, if the probability that the three-site chain is in the 
configuration (-) — 0) is denoted P12, then the transition 
rate to state P13 = (H ) (corresponding to the ejec- 
tion of a proton from the last site into the right reservoir) 
is f3e v ~ H ~ K and the steady state current J = /3£-P 
(where the sum £ runs over all configurations that con- 
tain a proton at the last site), will contain the term 

p V-H-Kp i2 

Monte-Carlo simulations were implemented for rela- 
tively small (N = 10) systems by randomly choosing a 
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site, and making an allowed transition with the prob- 
ability £eKp(Ei - E f )/r max , where r max is the maxi- 
mum possible transition rate of the entire system. In 
the next time step, a particle is again chosen at ran- 
dom and its possible moves are evaluated. The cur- 
rents were computed after the system reached steady- 
state by counting the net transfer of protons across all 
interfaces (which separate adjacent sites and the reser- 
voirs) and dividing by N + 1. Physical values of J 
are recovered by multiplying by r max . Particle occu- 
pation statistics within the chain were tracked by using 
the definitions of +,0 and — particle densities at each 
site i: p+(i) = (a^ai + l)/2), p {i) = ((1 - of)), and 
P-(i) — (cri{<Ji — l)/2), respectively. However, for our 
subsequent discussion, it will suffice to analyze simply the 
chain-averaged proton concentration ao = X)t=iPo(0- 
All MC results were checked and compared with the ex- 
act numerical results from the three-site, 27-state master 
equation. 

RESULTS AND DISCUSSION 

Here, we present MC simulation results for a lattice 
of size N = 10. The mechanisms responsible for the dif- 
ferent qualitative behaviors are revealed and the effects 
of each interaction term will be systematically analyzed. 
We explore a range of relative kinetic rates, all nondimen- 
sionalized in units of po, the intrinsic proton hopping rate 
from between adjacent waters. Estimates for po derived 
from quantum MP simul at ion s are on the order of lps -1 
ISadgeehi fc Cheng 19991 iMavri fc Berendsen 1991 
iMei et al. 1998llSchmitt fc Voth 1999j . 

One of the main features we wish to explore is the 
effect of multiple proton occupancy on current-voltage 
relationships. To understand what values of transition 
rates would permit multiple proton occupancy, consider 
water at pH=7, which has 10 _7 M proton and hydrox- 
yls. This concentration corresponds to about 60 H 3 + 
and 60 OH~ species per cubic micron. Even at pH 4, 
one would only have ~ 60, 000 hydroniums per /im 3 , cor- 
responding to a typical distance between hydroniums of 
~ 25nm. Since there are only ~ 10 — 20 waters within 
a single-file channel, and at pH 4, only about one in 
500,000 waters are protonated in bulk, multiple protons 
in a single channel can occur only if protonated species 
within the channel are highly stabilized by interactions 
with the chemical subgroups comprising the pore inte- 
rior. This stabilizing effect is modeled by small escape 
rates /?o > 7o i and assumed to be distributed equally such 
that po remains constant across all sites within the lat- 
tice. Although from a concentration point of view, small 
entrance rates cyo,5q arise from infrequent protons that 
wander into the first site of the channel, their exit rates 
/3o,7o can be suppressed even more by their stabiliza- 
tion once inside the channel. In all of our simulations, 
we will assume proton stabilization is moderately strong 
and limit ourselves rates /3o>7o < ao,<5o- The values we 
use give steady-state proton occupancies across the whole 
range of values from < 1 to N. 
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FIG. 4: Saturation due to small flip rates fc+ = fc_ = ko. Cur- 
rents and rates in all plots are nondimensionalized by units of 
po. (A) Small fco determines the rate limiting step whereupon 
increasing V does little to increase the current. Increasing fco 
pushes the sublinear (saturation) regime of the J — V rela- 
tionship to larger values of voltage V . (B) The total proton 
occupancy decreases with decreasing fco. 



First consider symmetric solutions and featureless, uni- 
form pores where ao = ^OiPa = 70 ■ The only possible 
driving force is an external voltage V. In Fig. 0] we 
plot the current-voltage relationship for various flipping 
rates fco- We initially ignore interaction effects and set 
H = K = R = 0. Currents for sufficiently small V 
are always nearly linear. However, for sufficiently large 
V, the rate limiting step eventually becomes the water 
flipping rate fco- Further increases in V do not increase 
the overall steady-state current, and the current-voltage 
curve becomes sublinear before saturating. The crossover 
to sublinear (water flipping rate limited) behavior de- 
pends on the value of fco > with sublinear onset occurring 
at higher voltages V for larger fco. In the noninteracting 
case, for most reasonable values of rate constants, any 
possible superlinear regime does not arise as it is washed 
out by the sublinear, water flip rate-limited saturation. 
The only instance found where noticeable superlinear be- 
havior in the steady-state proton current arises is in the 
limit of large fco and when ao , So , po <C /?o > 7o • For the pa- 
rameters explored, the currents J increase with increas- 
ing fco (Fig. 0J4); thus, the mean proton occupancies are 
qualitatively consistent with dynamics limited by internal 
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FIG. 5: Currents (A) and averaged proton occupation (B) in 
the presence of a constant water dipole-aligning field H > 0. 
For larger V, the ^-independent H assumption used in this 
scenario will break down due to the orientation effects of V 
on the water dipoles. 



proton hops. For small nipping rates, successive entry of 
protons is slow, while exit is not affected. As ka is in- 
creased, the bottlenecks near the entrance are relieved to 
a greater degree than those near the exit, increasing the 
overall proton occupancy (cf. Fig. 0]B). 

Figure[5]displays the effects of a fixed, external, dipole- 
orienting field H ^ 0. All other interactions and fields, 
except the external driving voltage V, are turned off. The 
convention used in the energy Eq. ^ favors a "+" state 
for H > 0. This asymmetry leads to an asymmetry in 
the J — V relationship (Fig. EJ4)- After an initial proton 
has traversed the channel, flipping of the "— " waters left 
in its wake is suppressed for H > 0, thereby preventing 
further net proton movement. The persistent blockade 
induced by increasing H is evident in Fig. \5jf3 where the 
proton density decreases for increasing H . 

Although H is held fixed in Fig. \5\ physically, dipole 
alignment fields arise from external electric fields that 
couple to the permanent dipoles of water. Therefore, we 
expect that H — LhvV where Lhv represents the orien- 
tational polarizability of the water molecule. It has been 
conjectured that when Lhv is positive (defined as prefer- 
ring waters with lone pairs pointing to the left, or in the 
"+" state) , the current should increase superlinearly with 
V since waters ahead of any proton will be oriented prop- 



FIG. 6: (A) Negative differential resistance (NDR) for large 

Lhv,V. Although transitions such as ... h0 h . . . — » 

. . . — h0 + + . . . are accelerated, giving rise to a state where 
proton transport to the right is possible, NDR can arise be- 
cause transitions such as . . . — h0 + + ...—>... — I — + ... 
created an additional — particle and is disfavored. (B) The 
average proton occupation decreases as V for large Lhv- 



erly as to receive it. Figure shows the current-voltage 
relationship for various Lhv- Although for very small 
Lhv, the current does increase very slightly, it becomes 
severely sublinear for larger Lhv an d V. In fact, it can 
attain a negative differential resistance (NDR) similar to 
that found in Gunn diodes or other "negistor" devices. 
The physical origins of NDR in proton conduction arise 
from the energetic cost of producing a "— " state as a 
proton moves forward. Although the path ahead of the 
proton is biased to "+" states, the proton transfer step 
as defined in our model necessarily leaves behind a "— " 
particle. Thus, although the field H = LhvV properly 
aligns waters ahead of a proton, it also provides an energy 
cost for the tail of "— " particles left by a forward-moving 
proton. This energetic penalty inhibits the proton from 
moving forward despite the direct driving force V acting 
on it. 

The average density plotted in Fig. decreases as V 
or large Lhv- Large Lhv not only hinders forward pro- 
ton hops, but enhances backward hops of protons that 
have just hopped forward during its previous time step. 
Proton dynamics are slowed dramatically, and only at 
the last site can they exit the pore. Proton entry from 
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a =8 =0.4; |3 ^y =0.05, k Q =2.0; H=K=0 
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FIG. 7: The effects of increasing nearest-neighbor proton- 
proton repulsion within the chain. Fixed parameters are ao = 
So = 0.4, (3 = 70 = 0.05, k = 2.0, and H = K = 0. (A) 
The onset of sublinear behavior in the J — V relationship 
is delayed for larger repulsions R, making the curves appear 
locally more superlinear. (B) The average proton densities 
per site. For small R, although densities are high, increasing 
V increases the clearance rate near the entrance such that the 
effectively increased injection increases overall proton density. 
At higher repulsions R, the clearance effects is not as strong 
and the simultaneously increased extraction rate prevents a 
large increase in the overall proton density. 

the left reservoir on the other hand, is often quickly fol- 
lowed by exit back into the left reservoir. The protons 
are effectively entry-limited, and the density is rather low. 
As V increases, the dynamics become even more "entry- 
limited," and the overall proton occupancy decreases. 

The effects of proton-proton repulsion (R > 0) are con- 
sidered in Figs and These simulations are consistent 
with the hypothesis that proto n-proton repulsions can 
give rise to superlinear current [Hille &: Schwarz 1978fl . 
Figure [7J4 shows a slight preference for superlinear be- 
havior as repulsion R is increased. Not surprisingly, Fig. 
0B shows that the overall density of protons within the 
pore decreases with increasing repulsion. 

The sublincar-to-supcrlincar behavior as the pro- 
ton concentration in the identical reservoirs is in- 
creased is shown in Fig. |HJ4. Although for these 
parameters, the effect is not striking, there is in- 
deed a trend away from sublinear behavior as pH 
is decreased, or, as ao — So is increased. Mea- 



R=4.0; R=y =0.05, k =2.0 




FIG. 8: Transition from sublinear to superlinear current be- 
havior as proton concentration in the symmetric reservoirs is 
increased. (A) J — V relationship for various concentrations 
a = 8 for fixed H = K = 0, R = 4.0, Po = 7o = 0.05, and 
ko = 2.0. (B) The averaged proton concentration ao at each 
lattice site as a function of driving voltage. The concentra- 
tions increase for all ranges of V as a = S is increased. 



surements, t hough, also show rather modest superlin- 
ear behavior lEisenman et al. 1980l iPhillips et al. 19991 
iRokitskava et al. 2002| ~ The occupancy also increases 
with decreasing pH, enhancing the effect of proton-proton 
repulsion. Thes e behaviors are consis tent with experi- 
mental findings [Eisenman et al. 1980j and those in the 
simulations depicted in Fig. where increased repulsion 
exhibited superlinear J — V curves. 

Finally, we consider the effects of dipole coupling K ^ 
between adjacent water molecules. This interaction is 
analogous to a nearest neighbor ferromagnetic coupling 
in e.g., Ising models. Fig. |5J4 shows that for sufficiently 
large ao — 5o, & superlinear behavior arises (for small 
enough V and large enough ko such that saturation has 
not yet occurred). Notice that as ao = do is increased, 
the J — V relationship can become more sublinear be- 
fore turning superlinear. Here, we have used a higher 
value of ko to suppress sublinear behavior to larger V, 
but the qualitative shift from sublinear to slightly super- 
linear behavior exists for small ko- Moreover, recent com- 
parisons between gramicidin A and gramicidin M chan- 
nels suggest that w ater reorientation is not rate-limiting 
|Gowen et al. 20 02], The nature of the superlinear be- 
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H=R=0; k=5; K=2.5; p =y =0.05 




FIG. 9: (A) The current-voltage relationship for various pro- 
ton injection rates in the presence of ferromagnetic water 
dipole coupling. (B) Mean proton occupations increase with 
increasing injection rates. 



havior can be deduced from Fig. [5]B, where the mean 
proton density is shown to increase with ota = Sq. Wa- 
ters that neighbor a proton are relieved of their dipolar 
coupling and can more readily flip to a configuration that 
would allow acceptance of another proton. For example, 
the transition ...0 — 0. + 0... will occur faster 

than ... ...—»■... — h0 . . .. This lubrication effect 

arises only when the proton density is high and K =/= 0. 

SUMMARY AND CONCLUSIONS 

We have developed a lattice model for proton conduc- 
tion that quantifies the kinetics among three approxi- 
mate states of the individual water molecules inside a 
simple, single-file channel such as gramicidin A. The 
three states represent water molecules with left and right- 
pointing water dipoles, and protonated ions. Our ap- 
proach allows us to explore the steady-state behavior of 
proton currents, occurring over timescales inaccessible 
by MD simulations. The model, along with analyses 
of Monte-Carlo simulations, also extends analytic mod- 
els [Schumaker et al. 200fj ISchumaker et al. 20 01] to in- 
clude multiple proton occupancy and the memory effects 
of protons that have recently traversed the water-wire. 
Monte-Carlo simulations of the lattice model was per- 
formed to test conjectures on a number of observed qual- 
itative features in proton transport across water wires. 



Four interaction energies that modify the kinetic rates 
are considered: A dipole-orienting field which tends to 
align the water molecules, a ferromagnetic dipole-dipole 
interaction terms between neighboring water molecules, a 
penalty from the repulsion between neighboring protons, 
and a external electric field (transmembrane potential) 
that biases the hops of the charged protons. 

We find current-voltage relationships that can be both 
superlinear and sublinear depending on the voltage V . 
For large enough voltages, the proton hopping step is no 
longer rate limiting. Water flipping rates limit proton 
transfer and further increases in V do little to increase 
the steady-state proton current J. This observation sug- 
gests that the observed transition from sublinear to su- 
perlinear behavior can be effected by varying an effective 
water flipping rate. Although we find that indeed proton- 
proton repulsion can lead to slightly superlinear J — V 
characteristics, particularly for large repulsions and pro- 
ton injection rates (low pH). 

Dipole-dipole interactions between neighboring wa- 
ters are also incorporated. Previous single-proton theo- 
ries iSchuma ker et al. 20001 [Schumaker et al. 200lj have 
considered the propagation of a single defect back and 
forth in the pore. In our model, the number of pro- 
tons and defects are dynamical variables that depend on 
the injection rates and the dipole-dipole coupling, respec- 
tively. For large coupling K, we expect very few defects, 
and effective water flipping rates will be low. However, 
when injection rates and proton occupancy in the pore is 
high, some dipole-dipole couplings are broken up by the 
intervening protons. Thus, protons can "lubricate" their 
neighboring dipoles, allowing them to flip faster than if 
they were neighboring a dipole pointed in the same direc- 
tion. Using simulations, we showed that this lubrication 
effect can give rise to a superlinear J — V relationship 

Although the parameters used in our analyses can 
be further refined by estimating them from shorter 
time MD simulation s , or other continuum app roaches 
[Edwards et al. 20021 iPartenskii fc Jordan 199^ . More 
complic ated local interactions with membrane lipid 
dipoles |Rokitskava et al. 2002j and internal pore con- 
stituents (such as Trp side groups |Dorigo et al. 19991 
iGowen et al. 20 02]) can be incorporated by allowing 
H, K,po and/or fco to reflect the local molecular environ- 
ment by var ying along the lat tice site (position) within 
the channel [Kolomeiskv 1998|] . 



The author thanks Mark Schumaker for vital discus- 
sions and comments on the manuscript. This work was 
performed with the support of the National Science Foun- 
dation through grant DMS-0206733, and the National 
Institutes of Health through grant R01 AI41935. 
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APPENDIX A: NOINTERACTING MEAN-FIELD 
RESULTS 

For the sake of completeness, and as a guide to aid 
qualitative understanding, we review analytic results in 
the case R = K = H = 0, where only exclusions are 
included. Some of these results have be en derived p revi- 
ously using mean-field approximations [Chou 2002j . 

If 7 = (£ = £o)) on ly pH differences between the two 
reservoirs can affect a nonzero steady-state proton cur- 
rent. The proton concentration difference is reflected by a 
difference between the entry rates from the two reservoirs 
Q'o 7^ <5oj aud the steady-state current can be expanded 
in powers of 1/N: J = ai / N + a 2 /N 2 + O (N~ 3 ) . In the 
long chain limit, we found |Chou 20 02] 



J 



In 



fc + fc_ 

N{k + +k-) A 

P(k + +k-)+k+8 7(fc + +fc-)+a(p-+fc-) 
7(fc_+fc + )+fc_a /3(fc + +fc_)+fc + <5(p_//c_+l) 



0(N- 2 ). 
(Af) 

For channels with reflection-symmetric molecular struc- 
tures, (3q = 70, and Eq. lAf I can be further simplified by 
expanding in powers of k-a — fc+<5, 



(3p + k-(k-a — k + S) 



N \J3(k- + k+) + k+5] (/? + S)(k+ + 

O ((k-a - k+S) 2 ) +0(1/N 2 ), 

Finally, in the large a and S — limit, 



(A2) 



J 



(k+ + k-)N 



log 1 



fc+ 
^k + k-P- 



aN(k+ + k-){k + + p_) 



(A3) 

For driven systems, where, say, a > <5, /3 > 7, and 
p + > p_, a finite current persists in the iV — » 00 
limit. We can use mean-field approximations famil- 
iar in the totally asymmetric simple exclusion process 
(TASEP) jDerrida f99llSchutz fc Domanv 1993] to con- 
jecture that three current regimes exist. If the both pro- 
ton entry and exit is fast, and the rate- limiting steps 
involve water flipping, or interior protons hops with rate 
p+, we expect that a maximal current regime exists and 
that the densities of the three states along the interior of 
a long chain are sp atially unifo rm. Mean-field analysis 
from previous work [Chou 20021 ] yields 



J 



2(p + fc_ — p-kj, 
iP++P-) 2 



(P- +P+) 



k- 



\Jk+ + k-y/k- + k + + p + + p- 



(A4) 



For a purely asymmetric process, p- — 0, and the current 
approaches the analogous maximal-current expression of 
the single species TASEP, 



J(p- = 0) 



P+k- 



4(fc_ 



O 



P+ 



(A5) 



except for the additional factor of fc_/(fc_ + k + ) repre- 
senting the approximate fraction of time sites ahead of a 
proton arc in the + configuration. These approximations 
neglect the influence of protons that have recently passed, 
temporarily biasing the water to be in a "— " configura- 
tion. Therefore, it is not surprising that these results are 
accurate only in the fc_ 3> p limit. 

A similar approach is taken when the currents are entry 
or exit limited. From the mean-field approximation of the 
steady-state equation for p± near the channel entry, 



P+Pap+ + k + p + - k-p- = 



dp- 
~dt 

d P+ 7 1 7 n 

— = -ap+ - k +P+ + k^p- = 0, 



(A6) 



where we have for simplicity set p- =7 = 0. Upon using 
normalization p- + po + p + = 1, and Eqs. IA6I we find 
the mean densities near the left boundary 



P+ = 



(a + k-)(p+ - a) 
p+{a + k- + k + ) ' 

k- (p — a) 
p+(a + k- + k + ) ' 



(A7) 



and the approximate entry rate-limited steady-state cur- 
rent 



J « P+P0P+ = ap + 



ak_ (1 — a/p) 
(a + k-+ k+) ' 



(A8) 



This result resembles the steady-state current of the 
low density phase in t he simple exclusion process 
[Derrida 199llChou 200^ . except for the factor fc_/(a + 
k- + fc+) representing the fraction of time the first site is 
in the + state, and able to accept a proton from the left 
reservoir. 

When the rate /3 is rate-limiting, we consider the mean- 
field equations near the exit of the channel 



dp- 

dt 

dp+ 
dt 



= f3p + k + p + - k-p- = 

= -P+P0P+ + k -P- - k +P+ = 0, 



(A9) 



and their solutions 
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/3(fc++ P+ -/3) p 
The exit-limited steady-state current is thus 



J w p> 



fc- + & 



/3(fc_ + fc+) 



(A10) 



(All) 



The results above are derived from mean-field assump- 
tions which neglect correlations in particle occupancy 
between neighboring sites. Although mean-field theory 



happens to give exact results for the simple exclusion 
process, the results above are only exact in the large 
k+/p + limit, as h as been shown by Monte-Carlo simula- 
tions | Chou 2002} . Only in this limit, where the memory 
of a previously passing proton is quickly er ased, are th e 
mean- field results quantitatively accurate |Chou 2002j |. 
Nonetheless, the mean-field calculations of the simplified 
system (H = K = R = 0) yields qualitatively correct 
results for the steady-state current, provides a connec- 
tion with well-known results of the TASEP, and gives 
an explicit qualitative description of the mechanisms at 
play. 



[Agmon 1995] Agmon, N. 1995. The Grotthuss mechanism. 
Chem. Phys. lett. 244:456-462. 

[Alberts et al. 1994] B. Alberts, D. Bray, J. Lewis, M. Raff, 
K. Roberts and J. D. Watson, Molecular biology of the cell, 
(Garland Publishing, New York, 1994). 

[Akeson & Deamer 1991] Akeson, M., and D. W. Deamer. 
1991. Proton conductance by the gramicidin water wire: 
Model for proton conductance in the FIFO ATPases? Bio- 
phys. J., 60:101-109. 

[Anderson 1983] Andersen, O. S. 1983. Ion movement 
through gramicidin A channels. Single-channel measure- 
ments at very high potentials. Biophys. J., 41:119-133. 

[Bala et al. 1994] Bala, P., B. Lesyng, and J. A. McCammon. 
1994. Applications of quantum-classical and quantum- 
stochastic molecular dynamics for proton transfer pro- 
cesses. Chem. Phys. 180:271-285. 

[Boyer 1997] Boyer, P. 1997. The ATP synthase - a splendid 
molecular machine. Annu. Rev. Biochem. 66:717-749. 

[Busath et al. 1998] Busath, D. D., C. D. Thulin, R. W. Hen- 
dershot, L. R. Phillips, P. Maughn, C. D. Cole, N. C. Bing- 
ham, S. Morrison, L. C. Baird, R. J. Hendershot, M. Cot- 
ten, and T. A. Cross. 1998. Non-contact dipole effects on 
channel permeation. I. Experiments with (5F-indole)Trp- 
13 gramicidin A channels. Biophys. J. 75:2830-2844. 

[Chernyshev & Cukierman 2002] Chernyshev, A. and S. 
Cukierman. 2002. Thermodynamic View of Activation En- 
ergies of Proton Transfer in Various Gramicidin A Chan- 
nels. Biophys. J. 82:182-192. 

[Chou 1998] Chou, T. 1998. How Fast do Fluids Squeeze 
Through Microscopic Single-File Channels? Phys. Rev. 
Lett. 80:85-89. 

[Chou 1999] Chou, T. 1999. Kinetics and thermodynamics 
across single-file pores: solute permeability and rectified 
osmosis. J. Chem. Phys. 110:606-615. 

[Chou & Lohse 1999] Chou, T. and D. Lohse. 1999. Entropy- 
driven pumping in zeolites and ion channels. Phys. Rev. 
Lett. 82:3552-3555. 

[Chou 2002] Chou, T. 2002. A spin flip model for one- 
dimensional water wire proton transport. J. Phys. A. 
35:4515-4526. 

[Chou 2003] Chou, T. 2003. Ribosome recycling, diffusion, 
and mRNA loop formation in translational regulation. 
Biophys. J. 85:755-773. 

[Cotten et al. 1999] Cotten, M., C. Tian, D. D. Busath, R. 
B. Shirts, and T. A. Cross. 1999. Modulating dipoles 



for structure-function correlations in the gramicidin A 
channnel. Biochemistry. 38:9185-9197. 
[Cukierman 1997] Cukierman, S., E. P. Quigley, and D. S. 
Crumrine. 1997. Proton conductance in gramicidin A and 
its dioxolane-linked dimer in different bilayers. Biophys. J. 
73:2489-2502. 

[Deamer 1987] Deamer, D. W. 1987. Proton permeation of 
lipid bilayers. J. Bioenerg. Biomembr. 19:457-479. 

[Decornez 1999] Dcornez, H., K. Drukker, and S. Hammes- 
Schiffer. 1999. Solvation and hydrogen- bonding effects on 
proton wires. J. Phys. Chem. A. 103:2891-2898. 

[Derrida 1998] Derrida, B. 1998. An exactly soluble non- 
equilibrium system: The asymmetric simple exclusion pro- 
cess. Physics Reports 301:65-83. 

[DeCoursey & Cherny 1994] DeCoursey, T. E., and V. V. 
Cherny. 1994. Voltage-activated hydrogen ion currents. J. 
Membr. Biol. 141:203-223. 

[Dorigo et al. 1999] Dorigo, A. E., D. G. Anderson, and D. D. 
Busath. 1999. Noncontact dipole effects on channel per- 
meation. II. Trp conformations and dipole potentials in 
gramicidin A. Biophys. J. 76:1897-1908. 

[Edwards et al. 2002] Edwards, S., B. Corry, S. Kuyucak, and 
S.-H. Chung. 2002. Continuum Electrostatics Fails to De- 
scribe Ion Permeation in the Gramicidin Channel. Bio- 
phys. J. 83: 1348 - 1360. 

[Eisenman et al. 1980] Eisenman, G., B. Enos, J. Hagghmd, 
and J. Sandblom. 1980. Gramicidin as an example of a 
single-filing ionic channel. Ann. N. Y. Acad. Sci. 339:8-20. 

[Fisher & Kolomeisky 1999] Fisher, M. E., and A. B. 
Kolomeisky. 1999. The force exerted by a molecular motor. 
Proc. Natl Acad. Sci. USA. 96:6597-6602. 

[Grabe & Oster 2001] Grabe, M. and G. Oster. 2001. Regu- 
lation of Organelle Acidity. J. Gen. Physiol. 117:329-343. 

[Gowen et al. 2002] Gowen, J. A., J. C. Markham, S. E. Mor- 
rison, T. A. Cross, D. A. Busath, E. J. Mapes, and M. F. 
Schumaker. 2002. The Role of Trp Side Chains in Tuning 
Single Proton Conduction through Gramicidin Channels. 
Biophys. J. 83:880-898. 

[Grotthuss 1806] Grotthuss, C. J. T. 1806. Sur la 
decomposition de l'eau et des corps qu'elle tient en 
dissolution a l'aide de l'electricite galvanique, Ann. Chim. 
58:54-74. 

[Hille & Schwarz 1978] Hille, B., and W. Schwarz. 1978. 
Potassium channels as multi-ion single-file pores. J. Gen. 
Physiol. 72:409-442. 



11 



[Hladky & Haydon 1972] Hladky, S. B., and D. A. Haydon. 
1972. Ion transfer across lipid membranes in the presence 
of gramicidin A. I. Studies of the unit conductance chan- 
nel. Biochim. Biophys. Acta. 274:294-312. 

[Hu & Cross 1995] Hu, W., and T. A. Cross. 1995. Trypto- 
phan hydrogen bonding and electrical dipole moments: 
functional roles in the gramicidin channel and implications 
for membrane proteins. Biochemistry. 34:14147-14155. 

[Hummer et al. 2001] Hummer, C, J. C. Rasaiah, and J. P. 
Noworyta. 2001. Water conduction through the hydropho- 
bic channel of a carbon nanotube. Nature 414:188-190. 

[Jordan 1984] Jordan, P. C. 1984. The total electrostatic po- 
tential in a gramicidin channel. J. Membr. Biol. 78:91-102. 

[Kalra et al. 2003] Kalra, A., S. Garde, and G. Hummer. 
2003. Osmotic water transport through carbon nanotube 
membranes. Proc. Natl. Acad. Sci. USA 100:10175-10180. 

[Karimipour 1999] Karimipour, V. 1999. Multispecies asym- 
metric simple exclusion process and its relation to traffic 
flow Phys. Rev. E 59:205-212. 

[Kolomeisky 1998] Kolomeisky, A.B. 1998. Asymmetric sim- 
ple exclusion model with local inhomogeneity. J. Phys. A: 
Math. Gen. 31:1153-1164. 

[Lanyi 1995] Lanyi, J. K. 1995. Bacteriorhodopsin as a model 
for proton pumps. Nature. 375:461-463. 

[Levitt et al. 1978] Levitt, D. G., S. R. Elias, and J. M. Haut- 
man. 1978. Number of water molecules coupled to the 
transport of sodium, potassium and hydrogen ions via 
gramicidin, nonactin or valinomycin. Biochim. Biophys. 
Acta. 512:436-451. 

[Lynden-Bell & Rasaiah 1996] Lynden-Bell, R. M., and J. C. 
Rasaiah. 1996. Mobility and solvation of ions in channels. 
J. Chem. Phys. 105:9266-9280. 

[MacDonald & Gibbs 1969] MacDonald, C. T. and J. H. 
Gibbs. 1969. Concerning the Kinetics of Polypeptide Syn- 
thesis on Polyribosomes. Bioploymers, 7:707-725. 

[Marx et al. 1999] Marx, D., M. E. Tuckerman, J. Hutter, 
and M. Parrinello. 1999. The nature of the hydrated excess 
proton in water. Nature. 397:601-604. 

[Mavri & Berendsen 1995] Mavri, J., and H. J. C. Berend- 
sen. 1995. Calculation of the proton transfer rate using 
density matrix evolution and molecular dynamics simu- 
lations: inclusion of the proton excited states. J. Phys. 
Chem. 99:12711-12717. 

[Mei et al. 1998] Mei, H. S., M. E. Tuckerman, D. E. Sagnell, 
and M. L. Klein. 1998. Quantum nuclear ab initio molec- 
ular dynamics study of water wires. J. Phys. Chem. B. 
102:10446-10458. 

[Nagle 1987] Nagle, J. F. 1987. Theory of passive proton con- 
ductance in lipid bilayers. J. Bioenerg. Biomembr. 19:413- 
426. 

[Nagle & Horowitz 1978] Nagle, J. F., and H. J. Morowitz. 
1978. Molecular mechanisms for proton transport in mem- 
branes. Proc. Natl. Acad. Set. USA. 75:298-302. 

[Nagle & Tristan-Nagle 1983] Nagle, J. F., and S. Tristram- 
Nagle. 1983. Hydrogen bonded chain mechanisms for pro- 
ton conduction and proton pumping. J. Membr. Biol. 74:1- 
14. 

[Partenskii & Jordan 1992] Partenskii, M. B., and P. C Jor- 
dan. 1992. Nonlinear dielectric behavior of water in trans- 



membrane ion channels: ion energy barriers and the chan- 
nel dielectric constant. J. Chem. Phys. 96:3906-3910. 

[Phillips et al. 1999] Phillips, L. R., C. D. Cole, R. J. Hen- 
dershot, M. Cotten, T. A. Cross, and D. D. Busath. 1999. 
Noncontact Dipole Effects on Channel Permeation. HI. 
Anomalous Proton Conductance Effects in Gramicidin. 
Biophys. J. 77:2492-2501. 

[Pomes & Roux 1996] Pomes, R., and B. Roux. 1996. Struc- 
ture and dynamics of a proton wire: A theoretical study 
of H+ translocation along the single- file water chain in the 
gramicidin A channel. Biophys. J. 71:19-39 

[Pomes & Roux 1998] Pomes, R., and B. Roux. 1998. Free 
energy profiles for H+ conduction along hydrogen-bonded 
chains of water molecules. Biophys. J. 75:33-40. 

[Prokop & Skala 1994] Prokop, P., and L. Skala. 1994. The- 
ory of proton transport along a hydrogen bond chain in 
an external field. Chem. Phys. Lett. 223:279-282. 

[Rokitskaya et al. 2002] Rokitskaya, T. I., E. A. Kotova, and 
Y. N. Antonenko. 2002. Membrane Dipole Potential Mod- 
ulates Proton Conductance through Gramicidin Channel: 
Movement of Negative Ionic Defects inside the Channel. 
Biophys. J. 82:865-873. 

[Sadgeghi & Cheng 1999] R. R. Sadcghi and H.-P. 
Cheng. 1999. The dynamics of proton transfer in a 
water chain. J. Chem. Phys. 111:2086-2094. 

[Sagnella et al. 1996] Sagnella, D. E., K. Laasonen, and M. L. 
Klein. 1996. Ab initio molecular dynamics study of proton 
transfer in a polyglycine analog of the ion channel grami- 
cidin A. Biophys. J. 71:1172-1178. 

[Sagnella & Voth 1996] Sagnella, D. E., and G. A. Voth. 
1996. Structure and dynamics of hydronium in the ion 
channel gramicidin A. Biophys. J. 70:2043-2051 

[Schemer 1985] Scheiner, S. 1985. Theoretical studies of pro- 
ton transfers. Acc. Chem. Res. 18:174-180. 

[Schmitt & Voth 1999] Schmitt, U. W., and G. A. Voth. 1999. 
The computer simulation of proton transport in water. J. 
Chem. Phys. 111:9361-9381. 

[Schreckenberg et al. 1995] Schreckenberg, M., A. Schad- 
schneider, K. Nagel, and N. Ito. 1995. Phys. Rev. E 
51:2939-2949. 

[Schumaker et al. 2000] Schumaker, M. F., R. Pomes, and B. 
Roux. 2000. A Combined Molecular Dynamics and Diffu- 
sion Model of Single-Proton Conduction through Grami- 
cidin. Biophys. J. 78:2840-2857. 

[Schumaker et al. 2001] Schumaker, M. F. and R. Pomes, and 
B. Roux. 2001. Framework Model For Single Proton Con- 
duction through Gramicidin, Biophys. J. 80:12-30. 

[Schutz & Domany 1993] Schiitz, G. and E. Domany. 1993. 
Phase Transitions in an Exactly Soluble One-Dimensional 
Exclusion Process. J. Stat. Phys. 72:277-296. 

[Syganow & von Kitzing 1999] A. Syganow and E. von Kitz- 
ing. 1999. (In)validity of the Constant Field and Constant 
Currents Assumptions in Theories of Ion Transport. Bio- 
phys. J. 76:768-781. 

[Wu & Voth 2003] Wu, Y., and G. A. Voth. 2003. A Com- 
puter Simulation Study of the Hydrated Proton in a Syn- 
thetic Proton Channel. Biophys. J. 85:864-875. 



Water alignment, dipolar interactions, and multiple proton occupancy during 

water-wire proton transport 



Tom Chou 

Dept. of Biomathematics and IPAM, Los Angeles, CA 90095-1766 
(Dated: September 22, 2003) 

A discrete multistate kinetic model for water-wire proton transport is constructed and analyzed 
using Monte-Carlo simulations. The model allows for each water molecule to be in one of three states: 
oxygen lone pairs pointing leftward, pointing rightward, or protonated (HsO + ). Specific rules for 
transitions among these states are defined as protons hop across successive water oxygens. We then 
extend the model to include water-channel interactions that preferentially align the water dipoles, 
nearest-neighbor dipolar coupling interactions, and coulombic repulsion. Extensive Monte-Carlo 
simulations were performed and the observed qualitative physical behaviors discussed. We find the 
parameters that allow the model to exhibit superlinear and sublinear current-voltage relationships 
and show why alignment fields, whether generated by interactions with the pore interior or by 
membrane potentials always decrease the proton current. The simulations also reveal a "lubrication" 
mechanism that suppresses water dipole interactions when the channel is multiply occupied by 
protons. This effect can account for an observed sublinear-to-superlinear transition in the current- 
voltage relationship. 
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INTRODUCTION 

The transport of protons in aqueous media and across 
membranes is a fundamental process in chemical re- 
actions, solvation, and pH regulation in cellular en- 
vironments [Alberts et al. 1994, Grabe & Oster 2001]. 
Proton transport in confined geometries is also rele- 
vant for ATP synthesis [Boyer 1997] and light trans- 
duction by bacteriorhodopsin [Lanyi 1995]. In this 
paper, we develop a lattice model for describ- 
ing proton transport in one- dimensional environ- 
ments. This study is motivated by numerous mea- 
surements of proton conduction across channels em- 
bedded in lipid membranes [Akeson & Deamer 1991, 
Busath et al. 1998, Cotten et al. 1999, Cukierman 1997, 
Deamer 1987, Eisenman et al. 1980]. Experiments are 
typically performed using membrane-spanning grami- 
cidin channels that are only a few Angstroms in diame- 
ter. This geometric constraint imposes a single-file struc- 
ture on the configurations of the interior water molecules 
[Hille & Schwarz 1978, Hladky & Haydon 1972]. 

Under the same electrochemical potential gradients, 
conduction of protons across ion channels occurs at a 
rate typically an order of magnitude higher than that of 
other small ions. This supports a "water- wire" mech- 
anism [Akeson & Deamer 1991, Nagle & Horowitz 1978, 
Nagle & Tristan-Nagle 1983, Nagle 1987], first proposed 
by Grotthuss [Agmon 1995, Grotthuss 1806]. Across a 
water- wire, protons are shuttled across lone pairs of water 
oxygens as they successively protonate the waters along 
the single-file chain. However, since the hydrogens are 
indistinguishable, any one of the hydrogens in a water 
cluster (e.g., any of the three hydrogens on a hydronium) 
can hop forward along the chain to protonate the next 
water molecule or cluster of water molecules (cf. Fig. 1). 
This mechanism naturally allows much faster overall con- 



duction of protons compared to other small ions which 
have to wait for the entire chain of water molecules ahead 
of it to fluctuate across the pore in order to traverse the 
channel. 

A peculiar feature of measured current-voltage re- 
lationships is a crossover from sublinear to super- 
linear behavior as the pH of the reservoirs is low- 
ered. Measurements by Eisenman [Eisenman et al. 1980] 
were carried out in symmetric solutions in the 1-3 pH 
range, and the results were recently reproduced by 
Busath et al. [Busath et al. 1998] and Rokitskaya et 
al. [Rokitskaya et al. 2002]. These experiments were 
performed using simple, relatively featureless gramicidin 
A (gA) channels. One leading hypothesis is that the 
nonlinear proton current- voltage relationships arise from 
the intrinsic proton dynamics within such simple chan- 
nels. Specifically, multiple proton occupancy and re- 
pulsion among protons within the channel may give 
rise to the observed nonlinearity [Hille & Schwarz 1978, 
Phillips et al. 1999, Schumaker et al. 2001]. 

There have been a number of recent theoretical 
studies of water-wire proton conduction. Extensive sim- 
ulations on the quantum dynamics of proton exchange 
in essentially small, representative water clusters in 
vacuum have been used to predict microscopic hop- 
ping rates between water clusters [Bala et al. 1994, 
Sadgeghi & Cheng 1999, Marx et al. 1999, 

Mavri & Berendsen 1995, Mei et al. 1998, 

Sagnella et al. 1996, Schmitt & Voth 1999]. Pomes and 
Roux [Pomes & Roux 1996] have performed classical 
molecular dynamics (MD) simulations on water-channel 
interactions, proton hopping, and water reorientation. 
They derive effective potentials of mean force describing 
the energy barriers encountered by a single proton within 
the pore. Since MD simulations are presently limited to 
only processes that occur over a few nanoseconds, none 
of these computational methods are efficient at probing 
very long time, steady-state transport behavior. On 
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a more macroscopic, phenomenological level, Sagnella 
et al. [Sagnella & Voth 1996] and Schumaker et al. 
[Schumaker et al. 2000, Schumaker et al. 2001] have 
considered a the long-time behavior of a single proton 
and dipole "defect" diffusing in a single-file channel. 
The parameters used in these studies, including effective 
energy profiles and kinetic rates, were derived from MD 
simulations. Although the basic underlying structure 
assumed by all of these transport models qualitatively 
resembles the Grotthuss mechanism, they have not 
addressed multiple proton occupancy. 

In this paper, we will explore the intrinsically non- 
linear proton dynamics along a single-file water-wire. 
We use a dynamical lattice model that defines the dis- 
crete structural states of the water-wire that approxi- 
mate the continuous molecular orientations. Although 
the lattice model provides a different approach from MD 
simulations, it is more amenable to analysis at longer 
time scales, yet is connected to the microphysics inher- 
ent in MD simulations provided a consistent correspon- 
dence between the parameters is made. Rather than 
enumerating all possible molecular configurations, our 
lattice approach is resembles that developed for molecu- 
lar motors [Fisher & Kolomeisky 1999], mRNA transla- 
tion [MacDonald & Gibbs 1969, Chou 2003], traffic flow 
[Karimipour 1999, Schreckenberg et al. 1995], and ion 
and water transport in single-file channels [Chou 1998, 
Chou 1999, Chou & Lohse 1999]. Here, the proton occu- 
pancy along the water- wire will be self-consistently deter- 
mined by the prescribed lattice dynamics. The parame- 
ters used in our model are transition rates among discrete 
states that in principle can be independently computed 
from relatively short-time MD simulations. Despite the 
approximations inherent in our discrete model, it quali- 
tatively treats the effects of proton-proton repulsion and 
water- water dipole interactions. 

MODEL AND METHODS 

Qualitatively, protons hop from oxygen to oxygen dur- 
ing transport. The successive hops clearly do not have 
to involve an individually tagged proton; in this respect, 
proton currents resemble electrical conduction in a con- 
ductor. Many measurements of proton conduction across 
membranes are performed on the gramicidin model sys- 
tem. The interior diameter of gramicidin A (gA) is ~ 3—4 
Angstroms and can only accommodate water in a single 
file chain. Although the number of water molecules in 
this chain is a fluctuating quantity, their dynamics in and 
out of the channel will be assumed to be much slower than 
that of their orientational rearrangements and proton 
hopping [Hummer et al. 2001, Kalra et al. 2003]. We 
can thus treat the water wire as containing a fixed, av- 
erage number of water molecules. Within typical trans- 
membrane channels, are N £s 8 — 26 single-file waters 
[Levitt et al. 1978, Wu & Voth 2003]. 

Figure 1A shows a schematic of our model. We 
first assume that each "site" along the pore is occu- 
pied by a single oxygen atom which may either be 



part of neutral water, H2O, or a hydronium (H3O" 1 ") 
ion. Although protonated oxygens in bulk are often 
associated with larger complexes such as H 5 Oj (Zun- 
del cation), or HgO^j" (Eigen cation), in confined ge- 
ometries, the formation of the larger complexes is sup- 
pressed [Lynden-Bell & Rasaiah 1996]. Furthermore, we 
will show that our discrete model depicted in Fig. 1A 
can incorporate the dynamics of reactive proton transfer 
among transient clusters by an appropriate redefinition 
of a lattice site to contain the entire cluster. 
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FIG. 1: (A) Schematic of an iV = 11, three-species exclusion 
model that captures the steps in a Grotthuss mechanism of 
proton transport along a water wire. For typical ion channels 
that span lipid membranes, iV ~ 10 — 20. The transition rates 
are labeled in (B) and in the legend. Water dipole kinks are 
denoted by thick lines. 



Neutral waters have permanent dipole moments and 
electron lone-pair orientations that can rotate thermally. 
For simplicity, we bin all water dipoles (hydrogens) that 
point towards the right as "+" particles, while those 
pointing more or less to the left are denoted "— " par- 
ticles. The singly protonated species H 3 + is hybridized 
to a nearly planar molecule. Therefore, we will assume 
that hydronium ions are symmetric with respect to trans- 
ferring a proton forward or backwards, provided the ad- 
jacent waters are in the proper orientation and there are 
no external driving forces (electric fields). Each lattice 
site can exist in only one of three states: 0, +, or — , corre- 
sponding to protonated, right, or left states, respectively. 
Labeling the occupancy configurations er, = { — 1,0, +1}, 
allows for fast integer computation in simulations. 
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In addition to proton exclusion, the transition rules are 
constrained by the orientation of the waters at each site 
and are defined in Fig. IB. A proton can enter the first 
site (i = 1) from the left reservoir and protonate the first 
water molecule with rate a only if the hydrogens of the 
first water are pointing to the right (such that its lone- 
pair electrons are left-pointing, ready to accept a proton) . 
Conversely, if a proton exits from the first site back into 
the left reservoir (with rate 7), it leaves the remaining 
hydrogens right-pointing. In the pore interiors, a pro- 
ton at site i can hop to the right (left) with rate p + (p-) 
only if the adjacent particle is a right(left)-pointing, un- 
protonated water molecule. If such a transition is made, 
the water molecule left at site i will be left (right) point- 
ing. Physically, as a proton moves to the right, it leaves 
a wake of — particles to its left. A left moving proton 
leaves a trail of + particles to its right. These trails of — 
or + particles are unable to accept another proton from 
the same direction. Protons can follow each other suc- 
cessively only if water molecules can reorient such that 
these trails of +'s or — 's are thermally washed out. Wa- 
ter reorientation rates are denoted fc± (cf. Figs. IB and 
2). Protons at the rightmost end of the water wire (at 
site i = TV), exit with rate ft, which is different from p + 
since the local microenvironment (e.g., typical distance 
to acceptor electrons) of the bulk waters that accept this 
last proton is different from that in the pore interior. 
From the right reservoir, protons can hop back into the 
water wire with rate S if a water in the "— " configura- 
tion is at site i = N. The entrance rates a and 6 are 
functions of at least the proton concentration in the re- 
spective reservoirs. Figure 2 shows a representative time 
series of the evolution of a specific configuration. The 
rate-limiting steps in steady-state proton transfer across 
biological water channels are thought to be associated 
with water flipping [Pomes & Roux 1998]. 

The lattice discretization for individual H3O" 1 " ions 
need not be interpreted literally. Larger complexes can 
be effectively modeled by reinterpreting p± , k± , and the 
basic unit of hopping for the proton. For example, if 
certain conditions obtain, where ions are predominantly 
two-oxygen clusters (H 5 Oj), we defined each pair of wa- 
ters as occupying a single lattice site, k± as an effec- 
tive reorientation time for the following pair of waters, 
and p± as the hopping rate to an adjacent oxygen lone- 
pair. The Grotthuss water-wire mechanism is qualita- 
tively preserved as long as the proper identification with 
the microphysics is made. 

All eight "parameters" used in our model (the 
rates p±, k±, a, ft, 7, S), can be related to measured 
bulk quantities or derived from short-time MD simu- 
lations. They are a minimal set and are equivalent 
to the numerous bulk parameters used in other mod- 
els [Schumaker et al. 2001], such as the bulk proton dif- 
fusion constant, water orientational diffusion constants, 
etc. Using similar MD approaches then, one should be 
able to approximately fix the parameters used in our 
model. For example, variations in the potential of mean 
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FIG. 2: A time series depicting a number of representative 
transitions obeying the dynamical constraints of our model. 
A proton (0) at site i can move to the right with rate p+ 
only if site i + 1 is occupied by a properly aligned (lone-pair 
electrons pointing to the left) water molecule (+). When 
a proton leaves site i to the right, it leaves behind a water 
in state "— ", with lone pair electrons pointing to the right. 
Protons at site i can also move to the left with rate p- if 
site i — 1 is a water in the "— " state. In this case, a water 
is left behind at behind site i in the "+" state. The neutral 
water molecules must flip (+ O — ) in order for a nonzero 
steady-state current to exist. 
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FIG. 3: (A — D) Energy differences between final and initial 
states which involve a change in ferroelectric coupling, net 
dipole moment, and repulsive interactions. (E) A representa- 
tive energy barrier profile for H = K = R = V = (dashed 
curves). The energy profile for H, K, R, V ^ for a transition 
between the states considered in (D) is shown by the thick 
solid curve. 
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force along the pore (resulting from interactions of the 
different species with the constituents of the pore in- 
terior) are embodied by site-dependent transition rates 
p± and k±. Thus, MD-derived potentials of mean force 
used in previous models can also be implemented within 
our lattice framework. Such effects of local inhomo- 
geneities in the hopping rates have been studied an- 
alytically and with MC simulations in related models 
[Kolomeisky 1998]. 

The basic model described above has been studied an- 
alytically in certain limits where exact asymptotic re- 
sults for the steady-state proton current J were de- 
rived [Chou 2002]. However, this study did not explic- 
itly include any interactions other than proton exclusion 
and proton transfer onto properly aligned water dipoles. 
Effects arising from forces such as repulsion between 
protons in close proximity, interactions between water 
dipoles and external electric fields, and dipolar coupling 
between neighboring waters need to be considered. 

In Fig 3A, a proton moves down the electric potential 
reducing the total enthalpy by V, and a right-pointing 
dipole is converted into a left-pointing dipole at an en- 
ergy cost of H. Since both initial and final states have 
adjacent, repelling protons, the repulsion energy R does 
not enter in the overall energy change. In Fig. 3P, a 
proton moves down the potential (— V), a "+" water is 
converted to a "— ", (+H), a dipole "domain wall" is re- 
moved (-K), and the repulsive energy between adjacent 
charged protons is relieved (-R). The representation of 
these nearest neighbor effects can be succinctly written 
in terms of the energy of a specific configuration 

N-l JV 

E[{oi}] = -K^2aiO i+1 - H^Oi+ 

i=l i=l 

^ E^(i - ^ 2 )(i - - v E^Ii i(i - 

(i) 

The H,K,R,V parameters used in P[{<Jj}] are all in 
units of fee T and represent 

• H: energy cost for orienting a water dipole against 
external field 

• K: energy cost for two oppositely oriented, adja- 
cent dipoles 

• R: repulsive Coulombic energy of two adjacent pro- 
tons 

• V: energy for moving a charged proton one lattice 
site against an external field. 

V as the change in potential that a proton incurs as it 
moves between adjacent waters. The total transmem- 
brane potential V mem i, rane = NV. The local dielectric 
environment across a channel can induce a spatially 
varying effective potential Vi<j<jv [Edwards et al. 2002, 
Jordan 1984, Partenskii & Jordan 1992, 

Syganow & von Kitzing 1999]. In this study, we 



neglect this variation and assume constant V across the 
lattice. 

In order to connect the quantities H, K, J?, and V to 
the rates a, (5, 7, &,p±, k±, we will assume the transitions 
occur over thermal barriers. Although barriers to pro- 
ton hopping may be small, we employ the Arrhcnius 
forms in order to obtain a simple relationship so that 
qualitative aspects of the effects of H, K, R, and V can 
be illustrated. Activation-energy-based treatments for 
conduction across gramicidin channels have been pre- 
viously studied [Chernyshev & Cukierman 2002]. When 
the more complicated interactions and external poten- 
tials are turned on, the effective transition rates £ = 
{a, /3, 7, S, p±, k±} on which we base our Metropolis 
Monte-Carlo become 



^oexp^— j, (2) 

where £0 = {ao, fio, 70, So,Po, ko} are rate prefactors 
when H,K,R,V and AE are zero. In defining Eq. 2, 
we have assumed that the energy barrier due to the 
difference AE = P[{er-}] - E[{oi}} ({cr^} and {erj are 
the final and initial state configurations, respectively) is 
evenly split between the barrier energies in the forward 
and backward directions. We use the convention that 
p+ = p- = po and k+ = k- = ko when V = and H = 0, 
respectively. The constraints and the state-dependent 
transition rates determined by Eqs. 1 and 2 completely 
define a nonequilibrium dynamical model which we study 
using MC simulations. Note that in the original model 
(Fig. 2) we do not assume transition barriers, but rather 
only that the dynamics are Markovian. 

We first gained insight into the dynamics by consider- 
ing numerical solutions to the full Master equation for a 
short three site (N = 3) channel. If we explicitly enu- 
merated all 27 = 3 3 states of the three site model, the 
Master equation for the 27 component state vector P is 

^ = MP(t), (3) 

where M is the transition matrix constructed from the 
rates £. In steady-state, the Pi are solved by inverting 
M with the constraint YlfLi Pi = 1- The steady-state 
currents are found from the appropriate elements in P, 
times the proper rate constants in the model. For exam- 
ple, if the probability that the three-site chain is in the 
configuration (-1 — 0) is denoted P12, then the transition 
rate to state Pi 3 = (H ) (corresponding to the ejec- 
tion of a proton from the last site into the right reservoir) 
is /3e v ~ H ~ K and the steady state current J = /3 Y^'i Pi 
(where the sum J] ■ runs over all configurations that con- 
tain a proton at the last site), will contain the term 
P V - H - K P 12 . 

Monte-Carlo simulations were implemented for rela- 
tively small (N = 10) systems by randomly choosing a 
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site, and making an allowed transition with the prob- 
ability £exp(£j - E f )/r max , where r max is the maxi- 
mum possible transition rate of the entire system. In 
the next time step, a particle is again chosen at ran- 
dom and its possible moves are evaluated. The cur- 
rents were computed after the system reached steady- 
state by counting the net transfer of protons across all 
interfaces (which separate adjacent sites and the reser- 
voirs) and dividing by N + 1. Physical values of J 
are recovered by multiplying by v max • Particle occu- 
pation statistics within the chain were tracked by using 
the definitions of +, and — particle densities at each 
site v. p+(i) = (<Ji(<Ji + l)/2},p (i) = ((1 - a})), and 
p~(i) = {(Ji{oi — l)/2), respectively. However, for our 
subsequent discussion, it will suffice to analyze simply the 
chain-averaged proton concentration ctq = X^iPo(*)- 
All MC results were checked and compared with the ex- 
act numerical results from the three-site, 27-state master 
equation. 

RESULTS AND DISCUSSION 

Here, we present MC simulation results for a lattice 
of size TV = 10. The mechanisms responsible for the dif- 
ferent qualitative behaviors are revealed and the effects 
of each interaction term will be systematically analyzed. 
We explore a range of relative kinetic rates, all nondimen- 
sionalized in units of po, the intrinsic proton hopping rate 
from between adjacent waters. Estimates for p derived 
from quantum MD simulations are on the order of lps -1 
[Sadgeghi & Cheng 1999, Mavri & Berendsen 1995, 
Mei et al. 1998, Schmitt & Voth 1999]. 

One of the main features we wish to explore is the 
effect of multiple proton occupancy on current-voltage 
relationships. To understand what values of transition 
rates would permit multiple proton occupancy, consider 
water at pH=7, which has 10 _7 M proton and hydrox- 
yls. This concentration corresponds to about 60 H 3 0+ 
and 60 OH species per cubic micron. Even at pH 4, 
one would only have ~ 60, 000 hydroniums per pm 3 , cor- 
responding to a typical distance between hydroniums of 
~ 25nm. Since there are only ~ 10 — 20 waters within 
a single-file channel, and at pH 4, only about one in 
500,000 waters are protonated in bulk, multiple protons 
in a single channel can occur only if protonated species 
within the channel are highly stabilized by interactions 
with the chemical subgroups comprising the pore inte- 
rior. This stabilizing effect is modeled by small escape 
rates Po,7o, and assumed to be distributed equally such 
that po remains constant across all sites within the lat- 
tice. Although from a concentration point of view, small 
entrance rates ao,^o arise from infrequent protons that 
wander into the first site of the channel, their exit rates 
Ab7o can be suppressed even more by their stabiliza- 
tion once inside the channel. In all of our simulations, 
we will assume proton stabilization is moderately strong 
and limit ourselves rates Po,7o < ciq,So- The values we 
use give steady-state proton occupancies across the whole 
range of values from < 1 to TV. 
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FIG. 4: Saturation due to small flip rates k+ = k- = ko- Cur- 
rents and rates in all plots are nondimensionalized by units of 
po- (A) Small ko determines the rate limiting step whereupon 
increasing V does little to increase the current. Increasing ko 
pushes the sublinear (saturation) regime of the J — V rela- 
tionship to larger values of voltage V. (B) The total proton 
occupancy decreases with decreasing ko- 



First consider symmetric solutions and featureless, uni- 
form pores where cto = <5o,A) = 7o- The only possible 
driving force is an external voltage V. In Fig. 4, we 
plot the current-voltage relationship for various flipping 
rates ko- We initially ignore interaction effects and set 
H = K = R = 0. Currents for sufficiently small V 
are always nearly linear. However, for sufficiently large 
V, the rate limiting step eventually becomes the water 
flipping rate ko- Further increases in V do not increase 
the overall steady-state current, and the current- voltage 
curve becomes sublinear before saturating. The crossover 
to sublinear (water flipping rate limited) behavior de- 
pends on the value of ko , with sublinear onset occurring 
at higher voltages V for larger ko- In the noninteracting 
case, for most reasonable values of rate constants, any 
possible superlinear regime does not arise as it is washed 
out by the sublinear, water flip rate-limited saturation. 
The only instance found where noticeable superlinear be- 
havior in the steady-state proton current arises is in the 
limit of large ko and when ao,So,Po <C /?o , 7o • For the pa- 
rameters explored, the currents J increase with increas- 
ing ko (Fig. AA); thus, the mean proton occupancies are 
qualitatively consistent with dynamics limited by internal 
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FIG. 5: Currents (A) and averaged proton occupation (B) in 
the presence of a constant water dipole-aligning field H > 0. 
For larger V, the ^-independent H assumption used in this 
scenario will break down due to the orientation effects of V 
on the water dipoles. 



proton hops. For small flipping rates, successive entry of 
protons is slow, while exit is not affected. As ko is in- 
creased, the bottlenecks near the entrance are relieved to 
a greater degree than those near the exit, increasing the 
overall proton occupancy (cf. Fig. 4B). 

Figure 5 displays the effects of a fixed, external, dipole- 
orienting field H ^ 0. All other interactions and fields, 
except the external driving voltage V, are turned off. The 
convention used in the energy Eq. 1 favors a "+" state 
for H > 0. This asymmetry leads to an asymmetry in 
the J — V relationship (Fig. 5 A). After an initial proton 
has traversed the channel, flipping of the "— " waters left 
in its wake is suppressed for H > 0, thereby preventing 
further net proton movement. The persistent blockade 
induced by increasing H is evident in Fig. 5B where the 
proton density decreases for increasing H. 

Although H is held fixed in Fig. 5, physically, dipole 
alignment fields arise from external electric fields that 
couple to the permanent dipoles of water. Therefore, we 
expect that H = LhvV where Lhv represents the orien- 
tational polarizability of the water molecule. It has been 
conjectured that when Lhv is positive (defined as prefer- 
ring waters with lone pairs pointing to the left, or in the 
"+" state), the current should increase superlinearly with 
V since waters ahead of any proton will be oriented prop- 



FIG. 6: (A) Negative differential resistance (NDR) for large 

Lhv,V. Although transitions such as ... h0 h . . . —¥ 

... — 1-0 + + ... are accelerated, giving rise to a state where 
proton transport to the right is possible, NDR can arise be- 
cause transitions such as . . . — h0 + + ... — > ... — | — + ... 
created an additional — particle and is disfavored. (B) The 
average proton occupation decreases as V for large Lhv- 



erly as to receive it. Figure 6 shows the current- voltage 
relationship for various Lhv- Although for very small 
Lhv, the current does increase very slightly, it becomes 
severely sublinear for larger Lhv and V. In fact, it can 
attain a negative differential resistance (NDR) similar to 
that found in Gunn diodes or other "negistor" devices. 
The physical origins of NDR in proton conduction arise 
from the energetic cost of producing a "— " state as a 
proton moves forward. Although the path ahead of the 
proton is biased to "+" states, the proton transfer step 
as defined in our model necessarily leaves behind a "— " 
particle. Thus, although the field H = LhvV properly 
aligns waters ahead of a proton, it also provides an energy 
cost for the tail of "— " particles left by a forward-moving 
proton. This energetic penalty inhibits the proton from 
moving forward despite the direct driving force V acting 
on it. 

The average density plotted in Fig. 6B decreases as V 
or large Lhv- Large Lhv not only hinders forward pro- 
ton hops, but enhances backward hops of protons that 
have just hopped forward during its previous time step. 
Proton dynamics are slowed dramatically, and only at 
the last site can they exit the pore. Proton entry from 
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FIG. 7: The effects of increasing nearest-neighbor proton- 
proton repulsion within the chain. Fixed parameters are ao = 
So = 0.4, /3 = 70 = 0.05, k = 2.0, and H = K = 0. (A) 
The onset of sublinear behavior in the J — V relationship 
is delayed for larger repulsions R, making the curves appear 
locally more superlinear. (B) The average proton densities 
per site. For small R, although densities are high, increasing 
V increases the clearance rate near the entrance such that the 
effectively increased injection increases overall proton density. 
At higher repulsions R, the clearance effects is not as strong 
and the simultaneously increased extraction rate prevents a 
large increase in the overall proton density. 

the left reservoir on the other hand, is often quickly fol- 
lowed by exit back into the left reservoir. The protons 
are effectively entry-limited, and the density is rather low. 
As V increases, the dynamics become even more "entry- 
limited," and the overall proton occupancy decreases. 

The effects of proton-proton repulsion (R > 0) are con- 
sidered in Figs 7 and 8. These simulations are consistent 
with the hypothesis that proton-proton repulsions can 
give rise to superlinear current [Hille & Schwarz 1978]. 
Figure 7^4 shows a slight preference for superlinear be- 
havior as repulsion R is increased. Not surprisingly, Fig. 
7B shows that the overall density of protons within the 
pore decreases with increasing repulsion. 

The sublinear-to-superlinear behavior as the pro- 
ton concentration in the identical reservoirs is in- 
creased is shown in Fig. 8A. Although for these 
parameters, the effect is not striking, there is in- 
deed a trend away from sublinear behavior as pH 
is decreased, or, as cvq = Sq is increased. Mea- 



R=4.0; p =7 =0.05, k =2.0 




FIG. 8: Transition from sublinear to superlinear current be- 
havior as proton concentration in the symmetric reservoirs is 
increased. (A) J — V relationship for various concentrations 
ao = <5o for fixed H = K = 0, R = 4.0, /3 = 70 = 0.05, and 
ko = 2.0. (B) The averaged proton concentration ao at each 
lattice site as a function of driving voltage. The concentra- 
tions increase for all ranges of V as a = S is increased. 



surements, though, also show rather modest superlin- 
ear behavior [Eisenman et al. 1980, Phillips et al. 1999, 
Rokitskaya et al. 2002]. The occupancy also increases 
with decreasing pH, enhancing the effect of proton-proton 
repulsion. These behaviors are consistent with experi- 
mental findings [Eisenman et al. 1980] and those in the 
simulations depicted in Fig. 7 where increased repulsion 
exhibited superlinear J — V curves. 

Finally, we consider the effects of dipole coupling K ^ 
between adjacent water molecules. This interaction is 
analogous to a nearest neighbor ferromagnetic coupling 
in e.g., Ising models. Fig. 9A shows that for sufficiently 
large ceo = So, a, superlinear behavior arises (for small 
enough V and large enough ko such that saturation has 
not yet occurred). Notice that as cvo = is increased, 
the J — V relationship can become more sublinear be- 
fore turning superlinear. Here, we have used a higher 
value of ko to suppress sublinear behavior to larger V, 
but the qualitative shift from sublinear to slightly super- 
linear behavior exists for small ko- Moreover, recent com- 
parisons between gramicidin A and gramicidin M chan- 
nels suggest that water reorientation is not rate-limiting 
[Gowen et al. 2002]. The nature of the superlinear be- 
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FIG. 9: (A) The current-voltage relationship for various pro- 
ton injection rates in the presence of ferromagnetic water 
dipole coupling. (B) Mean proton occupations increase with 
increasing injection rates. 



havior can be deduced from Fig. 9B, where the mean 
proton density is shown to increase with ao = Sq. Wa- 
ters that neighbor a proton are relieved of their dipolar 
coupling and can more readily flip to a configuration that 
would allow acceptance of another proton. For example, 
the transition ...0 — 0...— J-...0 + 0... will occur faster 

than ... . . . — > ... — HO .... This lubrication effect 

arises only when the proton density is high and K ^ 0. 

SUMMARY AND CONCLUSIONS 

We have developed a lattice model for proton conduc- 
tion that quantifies the kinetics among three approxi- 
mate states of the individual water molecules inside a 
simple, single-file channel such as gramicidin A. The 
three states represent water molecules with left and right- 
pointing water dipoles, and protonated ions. Our ap- 
proach allows us to explore the steady- state behavior of 
proton currents, occurring over timescales inaccessible 
by MD simulations. The model, along with analyses 
of Monte-Carlo simulations, also extends analytic mod- 
els [Schumaker et al. 2000, Schumaker et al. 2001] to in- 
clude multiple proton occupancy and the memory effects 
of protons that have recently traversed the water-wire. 
Monte-Carlo simulations of the lattice model was per- 
formed to test conjectures on a number of observed qual- 
itative features in proton transport across water wires. 



Four interaction energies that modify the kinetic rates 
are considered: A dipole-orienting field which tends to 
align the water molecules, a ferromagnetic dipole-dipole 
interaction terms between neighboring water molecules, a 
penalty from the repulsion between neighboring protons, 
and a external electric field (transmembrane potential) 
that biases the hops of the charged protons. 

We find current- voltage relationships that can be both 
super linear and sublinear depending on the voltage V. 
For large enough voltages, the proton hopping step is no 
longer rate limiting. Water flipping rates limit proton 
transfer and further increases in V do little to increase 
the steady-state proton current J. This observation sug- 
gests that the observed transition from sublinear to su- 
perlinear behavior can be effected by varying an effective 
water flipping rate. Although we find that indeed proton- 
proton repulsion can lead to slightly superlinear J — V 
characteristics, particularly for large repulsions and pro- 
ton injection rates (low pH). 

Dipole-dipole interactions between neighboring wa- 
ters are also incorporated. Previous single-proton theo- 
ries [Schumaker et al. 2000, Schumaker et al. 2001] have 
considered the propagation of a single defect back and 
forth in the pore. In our model, the number of pro- 
tons and defects are dynamical variables that depend on 
the injection rates and the dipole-dipole coupling, respec- 
tively. For large coupling K, we expect very few defects, 
and effective water flipping rates will be low. However, 
when injection rates and proton occupancy in the pore is 
high, some dipole-dipole couplings are broken up by the 
intervening protons. Thus, protons can "lubricate" their 
neighboring dipoles, allowing them to flip faster than if 
they were neighboring a dipole pointed in the same direc- 
tion. Using simulations, we showed that this lubrication 
effect can give rise to a superlinear J — V relationship 

Although the parameters used in our analyses can 
be further refined by estimating them from shorter 
time MD simulations, or other continuum approaches 
[Edwards et al. 2002, Partenskii & Jordan 1992]. More 
complicated local interactions with membrane lipid 
dipoles [Rokitskaya et al. 2002] and internal pore con- 
stituents (such as Trp side groups [Dorigo et al. 1999, 
Gowen et al. 2002]) can be incorporated by allowing 
H, K,pq and/or ko to reflect the local molecular environ- 
ment by varying along the lattice site (position) within 
the channel [Kolomeisky 1998]. 
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APPENDIX A: NOINTERACTING MEAN-FIELD 
RESULTS 

For the sake of completeness, and as a guide to aid 
qualitative understanding, we review analytic results in 
the case R = K = H = 0, where only exclusions are 
included. Some of these results have been derived previ- 
ously using mean-field approximations [Chou 2002]. 

If V = (£ = £0), only pH differences between the two 
reservoirs can affect a nonzero steady-state proton cur- 
rent. The proton concentration difference is reflected by a 
difference between the entry rates from the two reservoirs 
cxq ^ Sq, and the steady-state current can be expanded 
in powers of l/N: J = aJN + a 2 /N 2 + (iV~ 3 ). In the 
long chain limit, we found [Chou 2002] 



J 



In 



N(k + +k-) ' 



P(k + + k-)+k + S -y(k + +k- 



_y(k-+k+) + k-a f3(k++k-) + k+d(p- /k- + l) 

] +0(N-*). 

(Al) 

For channels with reflection-symmetric molecular struc- 
tures, f3o = 70, and Eq. Al can be further simplified by 
expanding in powers of — k + S, 



J 



f3p+k-(k-a — k + 6) 



- + 



N[8(k.+k+) + k+S](d + S)(k + + k.) (A2 ) 
O ((jfe_a - k+S) 2 ) + 0(1/N 2 ), 

Finally, in the large a and S = limit, 



J 



k+k- 



{k+ + fc_)iV 



log 1 + 



P- 



jk+k-P- 



aN(k + + k_){k + 



+ 0(a~ 2 N- 1 ). 

(A3) 

For driven systems, where, say, a > S, (3 > 7, and 
p+ > p-, a finite current persists in the TV — > 00 
limit. We can use mean-field approximations famil- 
iar in the totally asymmetric simple exclusion process 
(TASEP) [Derrida 1998, Schutz & Domany 1993] to con- 
jecture that three current regimes exist. If the both pro- 
ton entry and exit is fast, and the rate- limiting steps 
involve water flipping, or interior protons hops with rate 
p+, we expect that a maximal current regime exists and 
that the densities of the three states along the interior of 
a long chain are spatially uniform. Mean-field analysis 
from previous work [Chou 2002] yields 



2(p + k- -p-k + ) 
(P+ +P-) 2 



(p- +p+) 



+ k-+ k+- 



+ k- y/k- + k + + p + + p- 



(A4) 



For a purely asymmetric process, p_ = 0, and the current 
approaches the analogous maximal-current expression of 
the single species TASEP, 



J(P- = 0) 



P+k- 



4(fc_ + *+) 



+ 



E± 



(A5) 



except for the additional factor of k-/(k- + k + ) repre- 
senting the approximate fraction of time sites ahead of a 
proton are in the + configuration. These approximations 
neglect the influence of protons that have recently passed, 
temporarily biasing the water to be in a "— " configura- 
tion. Therefore, it is not surprising that these results are 
accurate only in the k + , fc_ ^> p limit. 

A similar approach is taken when the currents are entry 
or exit limited. From the mean-field approximation of the 
steady-state equation for p± near the channel entry, 



p+Pop+ + k+p+ - k-p- 



dp. 
dt 



(A6) 



where we have for simplicity set p- = 7 = 0. Upon using 
normalization p_ + po + p + = 1, and Eqs. A6, we find 
the mean densities near the left boundary 



_ (a + k-)(p+ - a) 
p+(a + k- + k + ) ' 



P+ 



k- (p — a) 
p + (a + k- + £4 



(A7) 



and the approximate entry rate-limited steady-state cur- 
rent 



J w p+poP+ = ap + 



ak- (1 — a/p) 
(a + fc_ +k+)' 



(A8) 



This result resembles the steady-state current of the 
low density phase in the simple exclusion process 
[Derrida 1998, Chou 2003], except for the factor k_/(a + 
k-+k + ) representing the fraction of time the first site is 
in the + state, and able to accept a proton from the left 
reservoir. 

When the rate (3 is rate-limiting, we consider the mean- 
field equations near the exit of the channel 



dp- 
~df 

dp+ 
dt 



Pp + k + p + - k-p- = 



-P+P0P+ + k-p- - k + p + = 0, 



(A9) 



and their solutions 
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p- 



13 

P+ = — • 

P+ 



The exit-limited steady-state current is thus 



Jvf3po 



k- + 



B(k. + k+) 



P+ 



(A10) 



(All) 



The results above are derived from mean-field assump- 
tions which neglect correlations in particle occupancy 
between neighboring sites. Although mean-field theory 



happens to give exact results for the simple exclusion 
process, the results above are only exact in the large 
k±/p± limit, as has been shown by Monte-Carlo simula- 
tions [Chou 2002]. Only in this limit, where the memory 
of a previously passing proton is quickly erased, are the 
mean-field results quantitatively accurate [Chou 2002]. 
Nonetheless, the mean-field calculations of the simplified 
system (H = K = R = 0) yields qualitatively correct 
results for the steady-state current, provides a connec- 
tion with well-known results of the TASEP, and gives 
an explicit qualitative description of the mechanisms at 
play. 
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